function res = computeMutualInformation(phonemes, phoneme_instances)


% incorporate ICA
% normalise variance of each dimension
% use histograms?


phoneme_ndata = zeros(1, numel(phonemes));
phoneme_pca_data = cell(1, numel(phonemes));
for p = 1:numel(phonemes)
    for i = 1:numel(phoneme_instances{p})
% % %         nvals = size(phoneme_instances{p}{i}.params, 1);
        phoneme_pca_data{p} = [phoneme_pca_data{p} phoneme_instances{p}{i}.params'];
    end
    phoneme_ndata(p) = size(phoneme_pca_data{p}, 2);
end

all_pca_data = [phoneme_pca_data{:}];
ndata = sum(phoneme_ndata);
mean_pca_data = mean(all_pca_data, 2);
vars_pca_data = var(all_pca_data, 1, 2);


% compute ICA decomposition
use_ica = 1;
if use_ica
    [A W] = fastica(all_pca_data, 'lastEig', 20, 'numOfIC', 20);
    all_ica_data = W*(all_pca_data - mean_pca_data*ones(1, ndata));
    phoneme_ica_data = cell(1, numel(phonemes));
    for p = 1:numel(phonemes)
        phoneme_ica_data{p} = W*(phoneme_pca_data{p}-mean_pca_data*ones(1, phoneme_ndata(p)));
    end
else
    A = [];
    W = [];
end


mean_ica_data = mean(all_ica_data, 2);
vars_ica_data = var(all_ica_data, 1, 2);


% % % % normalise data
% % % if 0
% % %     all_data = (all_data-vars_data*ones(1, ndata))./(sqrt(vars_data)*ones(1, ndata));
% % %     for p = 1:numel(phonemes)
% % %         phoneme_data{p} = (phoneme_data{p}-vars_data*ones(1, phoneme_ndata(p)))./(sqrt(vars_data)*ones(1, phoneme_ndata(p)));
% % %     end
% % % end
all_data = all_ica_data;
phoneme_data = phoneme_ica_data;


% compute mutual information
base_entropy = zeros(1, size(all_data, 1));
cond_entropy = zeros(1, size(all_data, 1));
for d = 1:size(all_data, 1)
    base_entropy(d) = gaussEntropy(mean(all_data(d, :)), var(all_data(d, :), 1));
    
    for p = 1:numel(phonemes)
        temp_prior = phoneme_ndata(p)/sum(phoneme_ndata);
        temp_entropy = gaussEntropy(mean(phoneme_data{p}(d, :)), var(phoneme_data{p}(d, :), 1));
        
        cond_entropy(d) = cond_entropy(d) + temp_prior*temp_entropy;
    end
end
[base_entropy; cond_entropy]
infos = base_entropy - cond_entropy;


res.phoneme_pca_data = phoneme_pca_data;
res.all_pca_data = all_pca_data;
res.mean_pca_data = mean_pca_data;
res.vars_pca_data = vars_pca_data;
res.phoneme_ica_data = phoneme_ica_data;
res.all_ica_data = all_ica_data;
res.mean_ica_data = mean_ica_data;
res.vars_ica_data = vars_ica_data;
res.A = A;
res.W = W;
res.base_entropy = base_entropy;
res.cond_entropy = cond_entropy;
res.mutual_infos = infos;



function e = gaussEntropy(mu, sigma2)

e = 0.5*(numel(mu) + numel(mu)*(2*pi) + log(det(sigma2)));
